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Abstract. 


This  paper  surveys  the  state  of  general  purpose  software  for  the 
solution  of  partial  differential  equations.  A discussion  of  the 
purported  capabilities  of  twenty-one  programs  is  presented.  No  testing 
of  the  routines  was  performed. 


1.  Introduction 


Mathematical  software  for  partial  differential  equations  (pde's)  is 
a relatively  new  area.  It  is  only  within  the  last  five  years  that  programs 
have  been  published  which  are  capable  of  solving  relatively  large  classes 
of  partial  differential  equations  using  effective  numerical  approximations. 

The  first  attempts  at  general  purpose  software  were  oriented  toward 
providing  programming  languages  that  facilitated  writing  the  equation  and 
its  approximation  (usually  finite  differences).  Two  examples  of  such 
languages  were  PDELAN  [10]  and  PDEL  [5]-  These  seem  not  to  have  been  too 
widely  used  and  are  currently  not  supported. 

Current  popular  software  now  is  produced  in  the  form  of  subroutires 
with  user-provided  subprograms  defining  the  equation.  There  are,  however, 
several  notable  exceptions  to  this  that  will  be  discussed  in  detail  below. 

This  paper  resulted  from  a reasonable  extensive  literature  search  under- 
taken to  determine  what  software  was  available.  The  guiding  principle  in  the 
selection  process  was  generality.  Hence,  no  coverage  has  been  provided  in 
the  very  large  special  application  areas  such  as  finite  element  codes  for 
structural  engineering,  nuclear  reactor  codes,  hydrodynamics  codes  (produced 
principally  at  government  weapons  laboratories ) , and  petroleum  reservoir 
modelling  codes.  Information  on  some  of  the  programs  in  these  areas  may  be 
obtained  from  the  Argonne  (National  Laboratory)  Code  Center  and  the  COSMIC 
program  library  at  the  University  of  Georgia  computation  center.  Also,  no 
attempt  has  been  made  to  determine  the  existence  of  codes  based  on  integral 
equation  techniques. 

Before  proceeding  further  this  caveat  must  be  tendered:  no  attempt  has 
been  made  to  ascertain  the  availability,  portability,  utility,  or  accuracy 
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of  any  of  the  software  mentioned  in  the  remainder  of  this  paper. 


2.  Mathematical  Preliminaries. 


In  this  section  are  discussed  some  mathematical  ideas  that  are  pertinent 
to  the  description  of  pde  software. 

For  time  dependent  pde's  the  most  popular  approximation  is  the  method  of 
lines  wherein  the  spatial  derivatives  are  approximated  in  some  manner  in  order 
to  yield  a system  of  ordinary  differential  equations  (.ode's).  The  time 
integration  of  this  system  of  ode's  is  then  accomplished  by  the  use  of  existing 
software  for  ode's. 

For  example,  suppose  we  wish  to  approximate  the  single  equation 

ut  = f (x,t,ux,uxx)  , 0 < x < 1 , t > 0 (2.1) 

with  boundary  conditions 

u(o,t)  = u (l , t ) =0  , t > 0 (2.2) 

and  initial  condition 

u(x,0)  = g(x)  , 0 < x < 1 . (2.3) 

We  may  choose  a finite  difference  grid  {x^ ) given  by 

x.  = ih  , i = 0,1,2. ..  ,N*-1  , h = -77  , 

1 7 n+1 

and  supposing  v^(t)  to  be  an  approximation  to  u(x^,t)  , we  replace  ux 

and  uxx  in  (2.1)  by,  say,  second-order  central  finite  difference  approxima- 
tions resulting  in  the  ode  system 

dv. 


V-t.  V . _ 

“dt  = f(xi't’(Tid'vi-i)/a,(vi-i-2VVi)/h  } 
= Fi(t'vi.i'vi'Vi)  7 


(2.4) 


for  i = 1,2 , . . . ,N  . From  (2.2)  we  have 

v0(t)  = Vnd(t)  = 0 

and  from  (2-3)  we  obtain  the  initial  conditions 

vi(0)  = g(xj  i = 1,2,. ..,N  . (2.5) 

Equations  (2.k)  and  (2.5)  form  an  initial  value  ode  problem  which  can  be 
integrated  by  some  ode  package. 

It  is  not  necessary  to  use  finite  difference  approximations  to  remove 
the  spatial  variation.  One  might  assume  that  the  solution  u of  (2.l)  may 
be  written 

n 

u(x,t ) = Z c-  (t )b.  (x) 
i=l  1 1 

where  the  known  functions  f b^(x)}  form  a suitable  basis.  Depending  on  the 
choice  of  method,  e.g.  collocation  or  Galerkin,  for  removing  the  spatial 
variation  one  obtains  an  initial  value  problem  for  a system  of  ordinary  dif- 
ferential equations  in  the  unknowns  {c.  (t)}  . 

Most  developers  of  pde  software  who  have  chosen  to  use  the  method  of  lines 
rely  heavily  on  existing  software  for  ode’s  and  for  the  computation  of  the 
basis  functions.  Specifically  the  ode  software  used  almost  exclusively  is  one 
of  the  several  codes  based  on  Gear’s  method  written  by  Hindmarsh  [ 12, 13, 1*4-]  • 

The  basis  functions  chosen  are  frequently  B-splines  and  the  package  of  de  Boor 
[2]  is  used. 

The  use  of  the  method  of  lines  for  two-space  dimension  parabolic  pde's 
introduces  a significant  complication  in  that  it  is  generally  accepted  that 
the  ode  system  to  be  solved  may  be  stiff.  Hence,  this  dictates  the  use  of 
implicit  methods  which  give  rise  to  very  large,  sparse  Jacobians  which  must  be 
inverted.  The  inversion  of  such  matrices  is  not  at  all  trivial  and  seems  to  be 
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the  principle  obstacle  in  the  development  of  multi-space  dimension  pde 

software.  This  same  problem  exists  for  the  development  of  robust  software 

for  elliptic  pde's  defined  over  general  regions. 

3-  Software  for  Hyperbolic  Systems 

Hyperbolic  equations  are,  by  far,  the  most  difficult  class  of  partial 

differential  equations  to  approximate  numerically.  The  principal  difficulties 
lie  in  determining  the  correct  approximation  at  the  boundaries  of  the  domain 
and  accurately  representing  discontinuities  in  the  solution. 

To  illustrate  the  former  difficulty,  suppose  we  want  to  approximate  the 
solution  of  the  simple  system 


Au 


O-U 


"u/ 

H 

1 

where  u = 

, A = 

f 

O 

» 

H 

0 < x < 1 


It  is  easily  seen  that  the  general  solution  is  u^(x,t)  = f(t*-x)  and 

Ug(x,t)  = g(t-x)  , where  f and  g are  arbitrary  functions  determined  by 
the  boundary  and  initial  conditions.  Clearly,  at  the  boundary  x = 0 , 

(if  t < l),  u^  must  be  specified  and  un  cannot  be  specified  unless  the 
solution  is  known.  At  x = 1 , u^  must  be  specified,  and  u^  cannot  be. 

However,  a straightforward  finite  difference  approximation  to  equation 
(3-l)  requires  that  two  conditions  be  given  at  each  boundary.  The  problem 
here,  of  course,  is  that  the  approximation  near  the  boundary  must  take  into 
account  the  nature  of  the  characteristics,  i.e.,  the  lines  x +•  t = constant 
and  x - t = constant.  At  x = 0 , the  characteristics  indicate  that  u^ 
should  be  computed  from  known  interior  values,  e.g.,  by  extrapolation  along 
the  outgoing  characteristics  x + t = constant.  Analogous  statements  apply 
at  the  boundary  x = 1 . 

Although  the  preceding  comments  may  seem  rather  obvious  for  the  system 
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(3-l),  where  A is  a constant  diagonal  matrix,  the  construction  of  accurate 
difference  approximations  becomes  more  difficult  when  A is  a function  of 
x,  t,  and  u . In  this  case  the  characteristics  must  be  computed  by  an 
eigenvalue/eigenvector  routine  at  each  time  step  and  at  each  boundary. 

3-1  DCG 

This  program,  written  by  Engquist  and  Smedsaas  [9]  , was  designed  to 
provide  software  for  the  numerical  solution  of  hyperbolic  systems  in  one 
space  dimension.  The  user  specifies  symbolically  the  equation 

ut  = f(x,t,u,ux)  , 

where  a < x < b , and  ueRn  . Initial  conditions,  boundary  conditions,  and 
certain  characteristics  of  the  problem,  e.g.  uniform  or  non-uniform  grid  in 
space,  or  that  some  particular  terms  should  be  considered  stiff  (and,  therefore, 
integrated  implicitly  in  time)  are  also'  specified.  These  specifications  are 
accepted  by  the  first  part  of  DCG,  the  analyser.  A syntactical  analysis  of  the 
equations  and  boundary  conditions  is  performed  to  detect  errors  in  the  specifica- 
tion and  the  problem  is  classified  by  linearity,  time  and  space  dependence  of 
coefficients,  etc.  The  principal  symbol,  i.e.,  the  matrix  J>f/fcux  , is 
examined  to  determine,  as  far  as  possible  symbolically,  whether  the  correct 
number  of  boundary  conditions  have  been  specified.  FORTRAN  code  is  prepared 
to  compute  the  eigenvalue  and  eigenvectors  during  execution,  if  necessary. 

(This  part  of  the  program  is  the  analysis  of  the  matrix  analogous  to  the 
matrix  A in  equation  (3.1)).  Difference  methods  and  boundary  strategies, 
i.e.  which  variables  should  be  extrapolated  at  the  boundaries,  are  selected 
according  to  the  preceding  analysis  and  user  specification,  and  a solution 
algorithm  is  determined.  This  algorithm  is  then  optimized. 

The  second  part  of  DCG,  the  synthesizer,  accepts  the  algorithm,  and, 
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using  a library  of  FORTRAN  codes,  constructs  a FORTRAN  program  to  solve  the 
given  problem.  This  program  contains  output  statements  that  save  all  inter- 
mediate execution  output  r - a disk.  After  execution,  a graphics  program 
operates  interactively  with  the  user  to  provide  graphs  of  any  of  the 
variables  as  functions  of  space  or  time. 

As  alluded  to  above,  the  approximation  is  done  by  finite  differences 
using  different  types  of  second  and  fourth  order  differencing.  Leapfrog, 
Crank-Nicolson,  or  semi-implicit  integration,  all  of  which  are  energy  con- 
serving. However,  the  user  may  specify  dissipation  if  desired. 

DCG  seems  to  be  a fairly  thorough  treatment  of  the  one- space -dimension 
problem,  but  its  availability  may  be  limited  by  the  fact  that  the  analyzer 
is  written  in  SIMULA  , a language  not  generally  available  in  numerical 
computation  centers. 


J.2  RKFPDE 

This  subroutine,  written  by  Gary  [11],  was  motivated  by  a desire  to 
provide  software  for  a class  of  problems  commonly  called  the  "stream  function - 
vorticity  equation" . Problems  in  this  class  occur  frequently  at  the  National 
Center  for  Atmospheric  Research,  where  this  software  was  developed. 

The  system  of  equations  is  assumed  to  have  the  form 


u,  = f + g + h(x,y,t,w,w  ,w  ,w  ,w  ,w  ,w  ,w  ) 
t x y x'  y xx'  xy  yy  xxxx  yyyy 


plus  an  optional  elliptic  equation 


Lp  = F(x,y,t,u,u  ,u  ,u  ,u  ,u  ) 
f > x ^ ^ xxxxJ  yyyy 


where  L is  assumed  to  be  separable,  and  where  u e Rn  , w = (u^jU^, • • • ,un,  p), 
f and  g are  functions  of  x,y,t,  and  w . The  domain  is  assumed  to  be  a 
rectangle  in  the  plane.  On  the  boundaries  of  the  rectangle  the  user  may 
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specify  one  of  the  conditions: 


1)  peri  .dicity 

2 ) symmetry  with  respect  to  the  boundary 

3)  skew  symmetry  with  respect  to  the  boundary 

h)  a mixed  condition,  or 

5)  extrapolation  of  the  characteristic  variables. 

(This  last  boundary  condition  requires  that  the  user  determine  those 
quantities  which  should  not  be  specified  at  a boundary). 

The  solution  is  accomplished  by  a finite  difference  approximation  on  a 
user-supplied  uniform  or  non-uniform  mesh.  The  method  of  lines  is  used.  The 
user  may  select  second  or  fourth  order  central  differences  in  space.  The  time 
integration  is  carried  out  by  a Runge-Kutta-Fehlberg  3-^  scheme  modified 
from  a routine  provided  by  Shampine.  The  elliptic  equation  is  solved  by  the 
subroutine  SEPELI,  a part  of  FISHPAK  (cf.  5*l)-  The  functions  f,g,h,  boundary 
data,  and  the  coefficients  of  the  elliptic  operator  are  given  in  user-supplied 
subprograms . 

This  software  was  written  for  the  Control  Data  7600  and  designed  to  use 
extended  core  memory.  However,  portability  should  not  be  too  affected  as  the 
calls  to  this  memory  have  been  isolated  in  the  code  and  may  be  replaced  by 
statements  appropriate  for  the  user's  computer.  A version  of  this  code  is 

currently  running  on  the  CRAY-1, 
h . Software  for  Parabolic  Systems. 

In  contrast  to  the  lack  of  general  purpose  codes  for  hyperbolic  systems 
there  are  a number  of  programs  available  for  parabolic  systems.  The  method 
of  lines  is  used  exclusively  in  these  codes.  Hence,  the  discussion  of  the 
approximation  in  each  code  will  be  separated  into  two  parts:  space  and  time. 
Tine  codes  will  be  discussed.  They  will  be  presented  in  an  order  corresponding 
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to  increasing  capabilities  of  the  programs.  (One  program,  GENEPI,  may 
also  be  used  for  solving  a parabolic  equation,  see  5-9) • 

4.1  DSs/g 

This  package,  produced  by  Schiesser  [28],  is  a successor  to  an  earlier 
code  hEANS  [27].  DSS/2  is  primarily  an  educational  tool.  It  was  designed 
as  a main  program  which  the  user  supplies  with  three  subroutines  that  define 
the  actual  pde  to  be  solved.  It  is  assumed  that  the  pde  can  be  written  in 
the  form 

u = F(x,t,u,u  ,u  ) , a<x<b,t>0  (4.1.l) 

U X XX 

where  u e Rn.  (The  expressions  u and  u mean  the  vectors  of  first  and 

x xx 

second  partial  derivatives,  resp. ) The  boundary  conditions  are  not  provided 
in  a separate  subroutine,  but  rather,  are  incorporated  by  the  user  into  the 
subroutine  that  defines  F . It  appears,  therefore,  that  only  relatively 

simple  boundary  conditions  may  be  handled  with  ease. 

In  order  to  evaluate  F , the  user  selects  from  a large  collection 

of  subprograms  within  DSS/2  the  one  that  automatically  computes  finite 

difference  approximations  to  u or  u . Centered  and  non-centered 

X XX 

differences  of  even  order  two  through  ten  are  available. 

The  time  integration  scheme  may  be  chosen  from  one  of  fourteen  different 
Runge-Kutta  one-step  methods  of  orders  one  through  five,  or  the  user  may 
select  the  Hindmarsh  GEARB  [14]  package. 

4.2  PDE PACK 

This  software  was  developed  by  Sincovec  and  Madsen  [22],  and  may  be 
purchased  or  leased  from  Scientific  Computing  Consulting  Services.  The 
package  consists  of  eleven  subroutines  with  communication  of  the  problem 
through  user-provided  subroutines.  The  system  of  pde's  is  assumed  to  be 
in  the  form 


8 


“ t'-  f<k)  (t'x’u'V  U°  \,1  ux1,J  *’ 


■K  '*\,»"xn)ix  1 ft-2-1’ 


where 

k = 1,2, . . . ,n,  a < x < b,  t > t , D = D (t,x,u), 

O X,J  K,J 

and.  c may  be  0,1,  or  2 . Terms  of  the  form 


- [ xC  D u ^ * ] 

c k,.i  x x 

X 

allow  the  program  to  handle  discontinuous  coefficients,  e.g.,  when  a problem  has 
material  interfaces.  The  inclusion  of  factors  xc  facilitates  the  treatment 
of  common  expressions  occurring  in  cartesian,  cylindrical,  or  spherical  coordinates. 

The  boundary  conditions  are  assumed  to  be  prescribed  at  each  boundary 
and  to  be  in  the  form 


. (k) 


(k) 

av  u'"'  + 3k  ux  = Yk  i k=l,2,...,  n 


If  ^ then  and  yk  may  te  functions  of  t and  u , but  if 

p = 0,  a k and  Yk  may  only  be  functions  of  t 

The  space  discretization  is  accomplished  via  central  finite  differences 
on  a grid  of  points  specified  by  the  user.  The  differences  are  second  order 
on  a uniform  grid,  but  only  first  order  on  a non-uniform  grid.  The 
differencing  has  been  designed  to  conserve  such  quantities  as  Du^  across 
material  interfaces.  The  finite-difference  grid  may  not  change  with  time. 

The  time  integration  is  accomplished  by  use  of  a modification  of 
the  Hindmarsh  package  GEAEB. 


k.3  PDECOL 

This  software  was  also  produced  by  Sincovec  and  Madsen  [20  ]•  The  user 
specifies  the  problem  via  a set  of  subroutines  and  calls  the  subroutine 
PDECOL  for  the  solution. 

The  equation  is  assumed  to  be  in  the  form 
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i 

u.  = f(t,x,u,u  ,u  ) , 

t X XX 

where  u c Rn  , a < x < b,  and  t > tQ  • The  boundary  conditions  are 
assumed  to  be  in  the  form 

•| 

b(u,u  ) = z(t)  . 

X 

However,  the  user  need  not  specify  a boundary  condition  at  all.  This  feature 
may  be  useful  if  the  user  attempts  to  solve  a hyperbolic  system  or  a coupled 
ode-pde  system. 

The  spatial  approximation  is  accomplished  by  the  use  of  collocation 
using  B-splines  as  the  basis.  The  breakpoints  for  the  B splines  are  provided 
by  the  user,  as  are  the  degrees  of  the  piecewise  polynomials  and  the  order 
of  continuity.  The  B-splines  are  computed  using  the  package  of  de  Boor  [2]. 

The  collocation  points  are  chosen  automatically  by  the  software. 

Since  the  resulting  system  of  ordinary  differential  equations  is  implicit, 
the  authors  use  GEARIB,  a modification  of  GEARS  developed  by  Hindmarsh  [lb]. 

k.4  HO  LID 

This  subroutine  package  was  developed  by  Hyman  [17].  As  in  the  previous 
two  codes,  the  user  writes  subroutines  defining  the  problem  and  then  calls 
T LID.  The  equation  must  be  in  the  form 

u,  = g(x,t,u,u  ,u  ,f  ) O+.U.l) 

t x xx  x 

where  u Rn  , a < x < b , t > t„  , and  f = f(x,t,u,u  ,u  ). 

0 x xx 

The  inclusion  of  the  argument  f on  the  right  side  of  (U.U.l'  allows 
conservative  differencing  of  advective  terms  which  ai-.se  naturally  in 
fluid  dynamics  equations. 

The  acceptable  boundary  conditions  are  given  by: 
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a)  periodicity  in  x, 

b)  au  ^ bu^  = c , where  a,b,  and  c may  be  functions  of  t and  u 
if  b f 0,  otherwise  a and  c may  only  be  functions  of  t , 

c)  some  differential  equation,  a free  boundary,  or  no  boundary 
condition  prescribed. 

The  spatial  discretization  is  accomplished  automatically  using  second  , 
fourth,  or  sixth  order  centered  differences,  fast  Fourier  transform 
approximation  (in  the  case  of  periodicity  only)  with  or  without  linear 
filtering  of  higher  modes,  or  unsymmetric  second  or  third  order 
differences . 

The  time  integration  is  carried  out  by  the  GEARB  package. 

Graphical  output  is  available.  The  plotting  procedures  use  standard 
CALCOMP  routines,  and,  as  such,  represent  a departure  from  the  portability 
standard  the  author  has  set. 

1+.5  ISQPDE 

This  software  by  Eason  and  Mote  [8]  is  designed  to  solve  the  equation 
ut  = F(x,t,u,u^,  u(2),  ...,  u^) 

n (i ) 

where  a < x < b , t > tQ,  u e R , and  u J represents  the  vector  of 
partial  derivatives  of  u of  order  j . At  the  boundaries,  general 
conditions  of  the  form 


G^  (x,t,u,u^,  . . . ,u^p_1^)  = 0 


x = a,b;  i = 1,2,. ..,q 


are  assumed  to  hold. 


The  approximation  technique  is  a combination  of  time  integration  and 
least  square  approximation.  The  user  selects  a set  of  points  in 





the  interior  and  on  the  boundary  of  [a,b]  at  which  points  the  eauation 

will  be  integrated  in  time.  A particular  set  of  basis  functions  (provided  by  the 

user)  is  used  to  represent  an  approximation  v(a.,x)  to  the  solution 

«3 

u(x,t.),  where  a.  represents  the  vector  of  coefficients  in  the  approximation 
J J 

at  time  t . . The  first  step  of  the  algorithm  is  to  form  a vector  of 
J 

residuals,  evaluated  at  the  points  [x^}  , for  initial  and  boundary  con- 
ditions. The  vector  aQ  , corresponding  to  the  initial  time  tQ  , is  determined 
by  minimizing  the  residual  in  the  least  squares  sense.  The  minimization  is 

accomplished  by  Powell's  non-linear  least  square  technique.  Assuming  a. 

J 

is  known,  the  approximation  v(a.,x)  is  used  in  a modified  Gear-Hindmarsh 

3 

pre  dictor-corrector  algorithm.  Firstly,  values  of  the  solution  are  predicted 
at  time  level  t.  ^ by  a Taylor  series  expansion  of  v(a.,x)  in  time. 

These  predicted  values  are  used  in  conjunction  with  the  corrector  equation 
to  form  a vector  of  residuals  at  the  points  [x^j  . This  vector  is  a 
function  of  the  unknown  coefficient  vector  a--*.-,  , which  is  determined  by  a 
least  squares  minimization  of  the  residual  vector. 

U.6  POST 

This  software,  written  by  Schryer  [29],  is  available  through  the  PORT 
library  of  Bell  Laboratories.  It  is  designed  to  solve  coupled  ordinary  -partial 
differential  equation  systems.  The  equation  treated  is 


f(x,t,u,ux,ut,uxt)  = a ( x,t,u,ut,u 


t xt  X 


where  a^x^b,  ue  R,  with  boundary  conditions 


b(t,u,ux,ut,uxt)  = 0 


at  x = a,b  . 


It  is  also  possible  to  impose  non-local  conditions  on  the  solution  such  as 
periodicity  or  that  the  integral  of  the  solution  over  some  interval  assumes 


i 


a given  value. 


The  spatial  variation  is  approximated  using  a Galerkin  technique  with  a 
B-spline  basis  (provided  by  the  de  Boor  package).  The  time  integration  is  a 
one-step  method,  explicit  or  implicit,  coupled  with  an  extrapolation  algorithm 
to  provide  automatic  error  control  in  the  time  integration. 


The  user  provides  a call  to  subroutine  POSTS  providing  information  about 
the  B-spline  basis,  time  integration  limits,  and  an  error  tolerance.  He  also 
provides  two  subroutines  for  the  evaluation  of  the  equation  and  boundary 
conditions. 

4.7  PDETWO 

PDETWO  was  written  by  Melgaard  and  Sincovec  [21].  It  is  a collection 
of  subroutines  that  provides  an  interface  between  a two  dimensional  pde  and 
an  ode  integration  package.  The  user  specifies  the  problem  by  providing  sub- 
routines and  calls  PDETWO  for  a solution. 

The  equation  is  assumed  to  have  the  form 


(k) 


u^ 


,(k) , 


(DHk,l 

u(l> 

X 

*x" 

u<n)) 

K,n  x x 

(DV,  , 

u(l) 

) 

..  , (DV  u(n) 

k,l 

y 

y 

k,n  y 

A 

A 

b2  ’ 

t > 

tQ  , k = 1,2, . . 

coefficients  DH  . and  DV  are  piecewise  continuous  functions  of 
k,j  k, J 

t,x,y  and  u . 

Boundary  conditions  are  assumed  to  have  the  form 


where  a^  , 
represents 


(k)  , (k)  . . _ 

a^u  +•  4>^un  = Ck  ’ k = l,2,...,n  , 

b^,  and  c^  are  piecewise  continuous  functions  and  u^ 

a derivative  normal  to  the  boundary. 
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The  user  provides  a two  dimensional  finite  difference  grid  and  the 


) 


i 


package  automatically  produces  an  approximation  at  each  grid  point  corresponding 
to  a five-point  star. 

The  time  integration  is  accomplished  by  the  use  of  GEARB.  The  use  of 
Newton's  method  to  solve  the  system  of  non-linear  equations  resulting  from 
the  use  of  implicit  integration  schemes  necessitates  the  solution  at  each 
Newton  iteration  of  a large  sparse  banded  system  of  equations.  The  coefficient 
matrix,  the  Jacobian  of  the  non-linear  system,  must  be  evaluated  at  each 
iteration  and  represents  a significant  computational  expense.  The  authors 
have  examined  various  techniques  for  this  task  and  have  in  their  view  implemented 
the  most  efficient  one. 


1+.8  DISPL 


This  software  package  is  the  collaborative  effort  of  Leaf,  Minkoff, 
Byrne,  Bleaknsy  and  Saltzman  [19]  of  the  Argonne  National  Laboratory.  The 
package  is  designed  specifically  for  one  and  two  spatially  dimensioned 
kinetics-diffusion  equations.  It  represents  an  effort  to  provide  reasonably 
general  purpose  software  designed  for  a particular  application. 

The  equation  is  assumed  to  thave  the  form 

[pCp]^(t,r,z,u)  u^k)  + 0 7 ' (vk(t,r,z,u)u('k^  )+  (l-©)vk(t,r,z,u)-V  u^ 

n n 


v< (D^(t,r,z,u,^u)  Vu 


ik 


c,  . u 

ki 


(i) 


dk: 


•1J 


i=l  i,j=l 


* ffe(t,r,z,u,7u)  , 

where  k = l,2,...,n  . The  equation  is  assumed  to  hold  over  a rectangle 
with  material  interfaces  in  the  (r,z) -plane.  To  increase  the  applicability 
of  the  code,  two  forms  of  advection  terms  (0=0  or  0 = l)  are  allowed. 

11+ 


The  boundary  conditions  on  the  sides  of  the  rectangle  are  assumed 
to  have  the  form 

, (k)  -*  „ (k)  -*  , .0 

a hu  + g ®kVu  ' n = y 

where  at  |3,  and  y may  depend  on  k and  the  side  of  the  rectangle,  n 
is  the  exterior  unit  normal,  and  p®  may  depend  on  u and  Vu  -|n|  . 

xC 

A Galerkin  procedure  is  used  as  the  spatial  approximation  using 
B-splines  as  the  basis  functions.  The  mesh  is  provided  by  the  user.  The 
B-splines  are  computed  using  de  Boor's  package. 

The  time  integration  uses  the  same  method  as  that  found  in  Hindmarsh's 
GEARIB  [13],  although  the  authors  do  not  use  GEARIB  itself. 

Extensive  graphics  capabilities  have  been  provided.  The  software  has 
been  written  in  IORTRAN  — a FORTRAN  preprocessor.  This  should  present  no 
problems  since  the  authors  are  including  in  the  external  distribution  a 
portable  copy  of  the  MORTRAN  translator. 

4.9  FORSIM  VI 

As  the  name  indicates,  the  FORSIM  package  of  Carver  et  al  [6]  is  in 
the  sixth  edition,  a testament  to  the  longevity  of  this  software  effort. 

This  package  contains  a main  program  plus  a large  number  of  subprograms 
that  provide  centered,  non-centered,  and  upwind  finite  difference  approximations 
of  various  orders  and  cubic  spline  approximations.  Time  integration  methods 
provided  are  a variable  step  Runge-Kutta-Fehlberg,  a fixed  step  fourth  order 
Runge-Kutta  with  Fehlberg  corrections,  a fixed  step  Euler,  and  the  Hindmarsh- 
Gear  algorithm.  For  this  last  algorithm  the  user  may  select  a number  of 
different  approximations  to  the  Jacobian,  ranging  from  no  Jacobian  (functional 
iteration)  to  the  full  Jacobian,  in  which  case  the  sparse  matrix  package  of 
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Curtis  -ind  Reid  [7]  is  used. 

To  use  this  package,  the  user  provides  a subroutine  UPDATE,  that 
contains  storage  declaration,  initial  conditions,  equation  and  boundary  condition 
specifications,  and  output  specifications.  To  define  the  pde  the  user  selects 
those  subprograms  corresponding  to  the  spatial  discretization,  creates  a 
section  of  code  which  performs  the  evaluation  of  the  right  side  of  the 
equation,  e.g.,  the  F in  (U.l.l),  and  incorporates  this  into  UPDATE.  The 
time  integration  is  performed  by  a call  within  UPDATE  to  another  subroutine. 

This  package  may  also  be  used  to  solve  two-  and  three-space  dimension 
problems  defined  on  rectangles  or  parallelepipeds,  respectively.  Graphical 
output  is  'also  available. 

5.  Joftware  for  Elliptic  Equations. 

Due  to  the  extensive  theory  that  exists  for  elliptic  partial  differential 
equations  and  the  large  amount  of  work  that  has  been  done  on  the  numerical 
approximation  of  the  solutions  of  such  equations,  it  is  not  surprising  to  see 
a large  selection  of  software  in  this  area.  However,  it  is  true  that  for 
complicated  problems,  e.g.  non-linear  or  non -separable,  the  solution  techniques 
are  not  as  reliable  and  software  is  not  generally  available. 

the  software  packages  will  be  described  in  an  order  corresponding  to 
increasing  complexity  of  the  type  of  problems  solved. 

5.1  FISHPAK  (Version  5) 

This  package  of  subroutines  was  written  by  Adams,  Swarztrauber  and 
Gweet  [l]  at  the  National  Center  for  Atmospheric  Research,  Boulder,  Colorado. 

It  is  a continuing  effort  brought  about  originally  in  response  to  the  need 
for  a complete  set  of  readily  available  well-tested,  reliable,  efficient,  and 
well-documented  software  to  solve  a subclass  of  elliptic  equations  which  occur 
frequently  in  the  study  of  geophysical  fluids. 
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The  subclass  consists  of  separable  elliptic  equations  with  particular 
emphasis  on  the  Poisson  equation  defined  on  a rectangle  in  cartesian, 
cylindrical  or  spherical  coordinates,  and  with  Dirichlet , Neumann  or 
periodic  boundary  conditions  prescribed.  Of  particular  importance  was  the 
need  to  automatically  treat  coordinate  singularities,  e.g.,  the  origin  in 
spherical  coordinates,  and  equation  singularities,  e.g.,  fully  periodic 
solutions. 

The  package  consists  of  eighteen  subroutines  and  one  subpackage. 

These  are: 

(a)  Twelve  drivers  that  define  second  order,  central  finite  difference 
approximations  on  staggered  and  unstaggered  uniform  grids , incorporate 
boundary  data,  and  treat  singularities  for  two-dimensional  modified 
Helmholtz  equations  in  cartesian,  polar,  cylindrical,  surface  spherical, 
and  spherical  cross-section  coordinate  systems,  for  a three-dimensional 
Helmholtz  equation,  and  a general  separable  two-dimensional  elliptic 
equation  without  coordinate  singularities, 

(b)  six  solvers  that  are  used  to  solve  the  linear  systems  of  equations 
arising  above, and  two  solvers  that  can  be  used  to  solve  finite  difference 
approximations  to  complex-valued  separable  elliptic  equations,  and 

(c)  a sub package  of  fast  Fourier  transform  routines  that  provide  periodic, 
sine,  cosine,  sine  quarter  wave,  and  cosine  quarter  wave  transforms  as  well 
as  the  full  complex  transform. 

The  solution  techniques  are  based  on  generalizations  of  the  Buneman 
variant  of  cyclic  odd-even  reduction  [3]  and  some  of  the  routines  provide, 
at  the  users  option,  fourth  order  approximations  using  the  method  of  deferred 
corrections . 

Buzbee  et  al  [k]  tested  version  1 of  the  package  and  made  criticisms 
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and  suggestions  [30]  to  the  authors.  These  suggestions  were  implemented  by 
the  authors  resulting  in  the  second  version  of  the  package. 

5-2  A Package  for  the  Helmholtz  Equation  on  an  Arbitrary  Region. 

This  package  of  four  subroutines  was  written  by  Proskurowski  [23]  at 
the  Lawrence  Berkeley Laboratories , Berkeley,  California.  It  may  be  used 
to  solve  the  Helmholtz  equation 

V2u  >■  cu  = f (5.1) 

defined  on  a general  bounded  planar  region  £1  with  either  Dirichlet  or 
Neumann  conditions  prescribed  on  the  boundary  of  the  region. 

The  package  solves  the  standard,  second  order,  central  finite  difference 
approximation  to  equation  (5-l)  by  embedding  the  region  in  a rectangle  R 
and  using  the  capacitance  matrix  technique  coupled  with  fast  Poisson  solvers. 

Briefly,  the  method  consists  of  four  steps: 

1)  generate  the  capacitance  matrix  C the  order  of  which  is  equal  to 
the  number  of  grid  points  within  fi  and  adjacent  to  the  boundary, 

2)  solve  (5.l)  on  R with  f arbitrarily  extended  to  R using  a fast 
Poisson  solver,  e.g.,  a routine  from  FISHPAK, 

3)  solve  the  capacitance  matrix  equation 

Cz  = b (5.2) 

and  use  z to  correct  the  values  of  f , and 

N)  solve  (5-l)  on  R with  the  corrected  f . 

For  a general  region,  step  3)  can  be  very  time  consuming  since  the  order 
of  C may  be  large.  This  also  may  make  the  direct  internal  storage  of 
C prohibitive. 
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To  efficiently  overcome  these  difficulties,  the  package  provides 
four  subroutines: 

1)  HLMHLZ  solves  equation  (5-2)  implicitly  via  a conjugant  gradient 
iteration,  thus  obviating  the  storage  of  C , 

2)  HELMIT  generates  and  stores  C explicitly  and  solves  (5-2) 
directly, 

3)  HELMSIX  solves  the  Dirichlet  problem  only  storing  C explicitly 
and  has  the  option  of  obtaining  fourth  or  sixth  order  approximations  via 
deferred  corrections,  and 

4-)  HELSYM  produces  a symmetric  approximation  to  the  Dirichlet  problem, 

using  an  explicitly  generated  C . This  routine  may  be  used  in  conjunction  with 

an  algorithm  for  the  computation  of  eigenvalues  and  eigenvectors  of  large  symmetric 

matrices,  e.g.,  the  Lanczos  method,  in  order  to  approximate  the  eigenvalues  of  the 
Laplacian. 

To  describe  the  irregular  boundary,  the  user  must  supply  a subroutine, 

DOMAIN,  that  specifies  the  coordinates  of  the  irregular  grid  points,  i.e. 
those  points  on  the  finite  difference  grid  for  which  some  of  its  neighbors 
are  not  within  the  domain.  The  user  must  specify  the  distance  from  the 
grid  point  to  the  boundary  and  the  associated  boundary  values.  The  right 
side  of  equation  (5-l)  is  furnished  in  the  user-supplied  subroutine  CHARGE. 

5-3  EIGEN 

Ryder  and  Sanderson  [26]  have  developed  a package  for  approximating  the 
eigenvalues  of  Laplace's  equation.  Specifically,  the  program  finds  approximate 
solutions  to 

V2u  + (u2u  =0  (5.3) 

defined  on  a bounded,  simply  connected  domain  in  the  plane,  subject  to 
Dirichlet  or  Neumann  boundary  conditions  on  various  portions  of  the  boundary 
of  the  domain. 
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A major  difficulty  in  the  approximation  is  the  treatment  of 
singularities  in  the  solution  induced  hy  re-entrant  corners,  i.e.  an 
interior  angle  greater  than  n , on  the  boundary.  They  overcome  this 
difficulty  by  seeking  solutions  to  equation  (5-3)  of  the  form 


k N. 

r-"1  r 


u 


(r,g)  = ) \ |c.  . J (our.  )sin  a-  ■ 6-  + d.  . J ^air.  )cos  a @ . \ 

L.  ' aqn-  1 iJ  1 iJ  aq,.  i ij  W 

i=l  j=0 


(5-J+) 


ij 


where  k is  the  number  of  corners  on  the  boundary,  are  determined 

from  the  angle  of  the  i—  corner,  and  (i\,e^)  are  'the  coordinates  of  the 

point  (r,Q)  in  a polar  coordinate  system  centered  at  the  i—  corner. 

Since  each  term  of  equation  (5*4)  is  a solution  of  equation  (5-3),  the 

coefficients  c. . and  d. . are  determined  by  attempting  to  satisfy  the 

ij 

given  boundary  condition  in  a least  squares  technique.  It  is  known  also 
that  such  an  expansion  of  the  solution  produces  the  correct  asymptotic 
singularity  at  each  corner. 

The  procedure  used  is  to  minimize  a certain  error  functional  in  cu  by 
selecting  an  uj  , solving  for  the  c„  and  d , then  trying  to  correct 
cu  by  minimizing  the  residual  of  the  least  squares  solution.  Such  local  minima 
are  "strong  candidates"  for  eigenvalues  even  though  convergence  proofs  for 
this  technique  do  not  exist.  The  routine  EIGEN  uses  the  routine  LOCALM 
written  by  R.  Brent  for  finding  local  minima,  and  the  routine  DSVD  written 
by  P.  Businger,  for  finding  the  singular  value  decomposition  of  the  matrix. 


5.E  itpack/region 

This  software  [18]  is  a result  of  the  continuing  research  at  the 
Center  for  Numerical  Analysis  of  the  University  of  Texas  at  Austin-  It  is 
a package  of  routines  which  approximate  the  solution  of  the  linear,  self- 
adjoint  elliptic  equation 
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(au  ) «-  (bu  ) + cu  = f (5-5) 

xx  y y 

defined  on  a somewhat  general  region  in  the  plane,  with  Dirichlet 
boundary  conditions  assumed  on  the  boundary  of  the  region.  The  coefficients 
a,b,  and  c may  be  functions  of  x and  y . 

The  solution  is  approximated  using  second  order  central  finite 
differences  defined  on  a uniform  grid  with  equal  spacing  in  x and  y . 

The  boundary  of  the  region  over  which  the  equation  is  assumed  to  hold  may 
consist  of  horizontal  or  vertical  grid  lines  or  lines  of  slope  +_  1 
connecting  grid  points.  The  region  may  have  holes  so  long  as  the  boundaries 
of  the  holes  satisfy  the  above  conditions.  Such  restrictions  eliminate 
the  occurence  of  irregular  grid  points  near  boundaries. 

To  solve  the  sparse  linear  system  arising  from  the  approximation  to 
equation  (5-5),  the  user  may  select  one  of  six  iterative  algorithms: 

1.  Jacobi  iteration  with  Chebyshev  acceleration. 

2.  Compresssed  Jacobi  iteration  with  conjugate  gradient  acceleration. 

3-  Jacobi  iteration  on  a reduced  system  with  Chebyshev  acceleration. 

U.  Jacobi  iteration  on  a reduced  system  with  conjugant  gradient  acceleration. 
5-  Symmetric  successive  over relaxation  with  Chebyshev  acceleration. 

6.  Symmetric  successive  overrelaxation  with  conjugant  gradient  acceleration. 
The  selection  of  the  acceleration  parameter  and  the  stopping  criterion  are 
automatic . 

The  subroutine  REGION  defines  the  grid  after  the  user  has  supplied  a 
polynomial  parameterization  of  the  boundary. 

5-5  POTENT 

This  piece  of  software,  written  by  Thomas  [32],  was  produced  to  aid 
engineers  in  the  solution  of  problems  in  electrostatics  and  magnetostatics. 
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fhe  equation,  written  in  divergence  form,  is 


v(e(x,y)vu)  - (eux)x  t-  ( uy)y  = f (x,y)  (5 .6) 

that  is  assumed  to  hold  cn  a general,  bounded  planar  region.  Dirichlet 
or  mixed  boundary  conditions  are  assumed  to  hold  on  the  boundary. 

The  method  of  approximation  is  to  use  a finite  difference  approximation 
to  equation  (5-6)  and  then  solve  the  resulting  system  by  ADI,  point  SOR,  or 
line  SOR  with  the  acceleration  parameter  either  specified  by  the  user  or 
set  internally  (by  assuming  Dirichlet  boundary  conditions  for  Laplace's  equation 
defined  on  a bounding  rectangle). 

The  user  must  provide  a subprogram  which  declares  whether  a point  is 
inside,  on,  or  outside  the  boundary,  and  , in  the  latter  case,  its  distance 
to  the  boundary.  Once  the  boundary  has  been  determined  it  is  outputed 
graphically  for  the  user  to  check  for  errors. 

5.6  ELLPACK  77 

This  software  research  project  is  coordinated  by  Rice  [24].  It  is  a 
cooperative  effort  among  many  people  interested  in  the  development  and 
evaluation  of  algorithms  related  to  solving  the  large,  sparse  linear  systems 
of  equations  arising  from  approximations  to  linear  elliptic  equations.  The 
goal  of  the  project  is  to  facilitate  the  testing  of  algorithms  that  only 
deal  with  a portion  of  the  total  solution  process.  This  is  done  by  defining 
a fixed  number  of  subproblems  (modules),  and  defining  fixed  interfaces 
between  these  modules.  There  are  currently  four  modules:  equation  formation, 
equation  indexing,  equation  solution,  and  output.  Researchers  may  contribute 
software  designed  for  one  of  the  modules.  More  importantly,  they  may  select 
existing  software  from  the  other  modules  for  testing  their  own  software, 
thereby  relieving  themselves  of  the  burden  of  unnecessary  coding.  The  interested 
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reader  is  referred  to  the  above  reference  for  a detailed  discussion  of 
the  ELLPACK  project. 

5-7  ELLCOL 

The  first  piece  of  software  from  the  ELLPACK  project  has 
been  designed  by  Houstis  and  Rice  [ 16 ] . It  is  designed  to  solve  the  linear 
equat ion 


defined  on  a general  two-dimensional  region  with  the  boundary  condition 


au  +-  bu 

x y 


cu  = g 


assumed  on  the  boundary.  The  coefficients  a,  g, . . ,a,b,  and  c may  be 

functions  of  x and  y 

The  method  of  approximation  is  to  embed  the  region  in  a rectangular 
grid  and  then  use  collocation  at  four  Gaussian  points  within  each  rectangle, 
using  bicubic  piecewise  Hermite  polynomials  for  the  basis.  The  approximation 
is  also  required  to  interpolate  the  boundary  condition  at  "appropriate”  points. 
The  system  is  solved  using  profile  Gauss  elimination. 

To  use  this  software,  the  user  must  provide  information  to  the  ELLPACK 
input  and  output  modules.  For  input  the  user  provides: 

1.  subroutine  BCOORD  that  provides  a parametric  representation  of  the 
boundary  of  the  region, 

2.  functions  COEF  and  BCOEF  that  provide  the  equation  and  boundary 
condition  coefficients,  and 

J.  function  F that  provides  the  functions  f and  g 
The  user  may  select  various  levels  of  output  of  intermediate  results, 


approximations  to  ux, 


u , and 

y 


and  the  execution  time. 
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5.8  EPDEl 


This  main  program  was  written  "by  Hornsby  [15]  at  CERN.  It  is  designed 
to  approximate  the  solution  of 

au  + bu  - cu  *■  du  + eu  = f 
xx  yy  x y 

where  a,b,...,f  may  be  functions  of  x,y,  and  u . The  equation  is 

assumed  to  hold  over  a region  bounded  by  a simple  closed  curves  C^,C  , . . , 

where  C , ...,C  must  lie  within  C . Along  the  curves  C.  the  following 
2 K -L  J 

boundary  conditions  are  assumed  to  hold: 

1.  Dirichlet , or 

2.  if  C.  is  a line  which  coincides  with  the  finite  difference  grid, 
then  a mixed  condition 

pu  +■  au  +■  ru  = s 
x * y 

may  hold. 

A uniform  finite  difference  grid  is  defined,  and  at  each  grid  point 
a finite  difference  equation  is  developed  using  the  four  neighboring  points. 

The  resulting  system  of  equations  is  solved  by  point  SOR.  The  acceleration 
parameter  oj  is  estimated  using  the  method  of  Carre*.  The  convergence  criterion, 
based  on  estimates  of  the  spectral  radius  of  the  iteration  matrix  involves  a 
percentage  accuracy  specified  by  the  user.  If  the  equations  are  non-linear, 
functional  iteration  is  coupled  with  the  SOR  iteration. 

The  input  for  this  program  is  complicated.  Not  only  does  the  user 
provide  subroutines  to  evaluate  the  coefficients  of  the  finite  difference 
equations,  but  the  user  must  define  each  boundary  by  a sequence  of  grid  points 
supplied  on  cards. 

5 -9  GENEPI 

This  software,  written  by  Roux  et  al  [25],  is  designed  to  solve  a general 
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linear  or  non-linear  elliptic  or  parabolic  equation  on  a two-dimensional 
rectangle  with  Dirichlet,  Neumann,  or  mixed  boundary  conditions.  This  program 
accepts  the  problem  definition  in  symbolic  form  and  generates  a FORTRAN 
program  that  solves  the  approximate  equations. 

The  user  may  specify  a uniform  or  non-uniform  grid  in  space  and  time. 

The  program  generates  the  finite  difference  approximation  using  the  four 
nearest  neighbors  and  solves  the  equations  using  point  or  line  relaxation 
or  ADI,  but  leaves  the  acceleration  parameter  selection  to  the  user. 

This  program  requires  disk  files,  so  portability  may  be  a problem. 

5-10  ELIPTI 

This  software,  written  by  Taylor  and  Taylor  [31] , is  designed  to 
solve  a general  non-linear  elliptic  equation  defined  on  an  arbitrary  region 
in  the  plane,  and  assuming  Dirichlet  boundary  conditions.  The  solution  is 
obtained  by  approximating  the  steady-state  solution  of  a related  time- 
dependent  parabolic  equation. 

The  approximation  is  by  finite  differences  on  a uniform  rectangular 
grid.  The  time  integration  is  performed  in  such  a way  that  the  scheme  is 
equivalent  to  an  ADI  scheme.  Aitken-Shanks  acceleration  is  used  in  time  also. 


The  author  would  like  to  express  his  appreciation  to 
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